Qries eduard-martinez.github.io

Qries @emartigo

Qries eduard-martinez

Qries efmartinez@icesi.edu.co

[1.] Imágenes satelitales

1.1 Motivación

Night lights and economic growth

Asian Financial Crisis: Java, Indonesia

Tomado de: Measuring Economic Growth from Outer Space
Tomado de: Measuring Economic Growth from Outer Space

Genocide Event: Rwanda

Tomado de: Measuring Economic Growth from Outer Space
Tomado de: Measuring Economic Growth from Outer Space

Puede acceder al documento aquí: Measuring Economic Growth from Outer Space

Measuring the size and growth of cities using nighttime light

Puede acceder al documento aquí: Measuring the size and growth of cities using nighttime light

Housing Vacancy Rate (VHR) y Luces Nocturnas

Puede acceder al documento aquí: An estimation of housing vacancy rate using NPP-VIIRS night-time light data and OpenStreetMap data

CAUSES OF SPRAWL: A PORTRAIT FROM SPACE

Urban Land in Miami, FL

Tomado de: CAUSES OF SPRAWL: A PORTRAIT FROM SPACE
Tomado de: CAUSES OF SPRAWL: A PORTRAIT FROM SPACE

Puede acceder al documento aquí: CAUSES OF SPRAWL: A PORTRAIT FROM SPACE

1.2 ¿Datos raster?

¿Qué es un raster?

Resolución

Bandas de un raster

1.3 Fuentes:

1.3.1 Night lights

DMSPL (1992-2013)

Puede acceder a los datos anuales aquí 1992-2103

VIIRS (2012-Currently)

  • Puede acceder a la web del proyecto aquí: link

  • Puede acceder a los datos mensuales aquí: link

  • Puede acceder a los datos diarios aquí: link

Harmonized (1992-2020)

Puede acceder a los datos y al paper aquí: link

1.3.2 Land cover

  • Puede acceder a la web del proyecto aquí: link

  • Puede acceder a los datos aquí: link

  • Puede acceder a la documentación del proyecto aquí: link ## Dictionary

1.3.3 NASA

Puede acceder a datos de lluvias, temperatura, contaminación, vientos y otros en la página oficial del proyecto GES DISC de la NASA aquí: https://disc.gsfc.nasa.gov.

[2.] Manipular datos raster en R

0.0 Configuración inicial

Para replicar las secciones de esta clase, primero debe descargar el siguiente proyecto de R y abrir el archivo clase-05.Rproj.

Instalar/llamar las librerías de la clase

## Instalar/llamar las librerías de la clase
require(pacman) 
p_load(tidyverse,rio,skimr,viridis,osmdata,
       raster,stars, ## datos raster
       ggmap, ## get_stamenmap
       spatstat, ## analis de clusters
       sf, ## Leer/escribir/manipular datos espaciales
       leaflet) ## Visualizaciones dinámicas
## 
## The downloaded binary packages are in
##  /var/folders/gd/6_fhny_d3y3c9l263sjz73dw0000gn/T//Rtmp4tzXcv/downloaded_packages
## solucionar conflictos de funciones
select = dplyr::select

2.1 Leer un raster

methods(class = "stars")
##  [1] [                 [[<-              [<-               %in%             
##  [5] $<-               adrop             aggregate         aperm            
##  [9] as_tibble         as.data.frame     as.im             as.owin          
## [13] as.POSIXct        c                 coerce            contour          
## [17] cut               dim               dimnames          dimnames<-       
## [21] droplevels        expand_dimensions filter            hist             
## [25] image             initialize        is.na             Math             
## [29] merge             mutate            Ops               plot             
## [33] prcomp            predict           print             pull             
## [37] rename            replace_na        select            show             
## [41] slice             slotsFromS3       split             st_apply         
## [45] st_area           st_as_sf          st_as_sfc         st_as_stars      
## [49] st_bbox           st_coordinates    st_crop           st_crs           
## [53] st_crs<-          st_dimensions     st_dimensions<-   st_downsample    
## [57] st_extract        st_geometry       st_geotransform   st_geotransform<-
## [61] st_interpolate_aw st_intersects     st_join           st_mosaic        
## [65] st_normalize      st_redimension    st_rotate         st_sample        
## [69] st_set_bbox       st_transform      st_write          time             
## [73] transmute         write_stars      
## see '?methods' for accessing help and source code
## importar raster de luces: raster
luces_r = raster('input/raster/night_light_202301.tif')
luces_r
## class      : RasterLayer 
## dimensions : 2990, 2919, 8727810  (nrow, ncol, ncell)
## resolution : 0.004166667, 0.004166667  (x, y)
## extent     : -79.01042, -66.84792, 0.002082733, 12.46042  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : night_light_202301.tif 
## names      : night_light_202301 
## values     : -0.64, 2208.49  (min, max)
## importar raster de luces: stars
luces_s = read_stars("input/raster/night_light_202301.tif")
luces_s
## stars object with 2 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
##                         Min. 1st Qu. Median      Mean 3rd Qu. Max.  NA's
## night_light_202301.tif  0.07    0.16   0.18 0.1810952     0.2 0.53 97553
## dimension(s):
##   from   to offset     delta refsys point x/y
## x    1 2919 -79.01  0.004167 WGS 84 FALSE [x]
## y    1 2990  12.46 -0.004167 WGS 84 FALSE [y]

2.2 Atributos de un raster

## resolucion
0.00416667*110000 
## [1] 458.3337
## geometria
st_bbox(luces_s)
##          xmin          ymin          xmax          ymax 
## -79.010415858   0.002082733 -66.847915761  12.460416166
st_crs(luces_s)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
st_dimensions(luces_s)
##   from   to offset     delta refsys point x/y
## x    1 2919 -79.01  0.004167 WGS 84 FALSE [x]
## y    1 2990  12.46 -0.004167 WGS 84 FALSE [y]
## atributos
names(luces_s)
## [1] "night_light_202301.tif"
names(luces_s) = "date_202301"

## valores del raster
luces_s[[1]] %>% max(na.rm = T)
## [1] 2208.49
luces_s[[1]] %>% min(na.rm = T)
## [1] -0.64
luces_s[[1]] %>% as.vector() %>% summary() 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      -1       0       0       0       0    2208 4031249
luces_s[[1]][is.na(luces_s[[1]])==T] %>% head()# Reemplazar NA's
## [1] NA NA NA NA NA NA
luces_s[[1]][2000:2010,2000:2010]
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24  0.28  0.14
##  [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16  0.25  0.14
##  [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10  0.17  0.12
##  [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12  0.14  0.12
##  [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15  0.13  0.11
##  [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15  0.10  0.13
##  [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14  0.13  0.13
##  [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13  0.19  0.19
##  [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09  0.13  0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11  0.11  0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18  0.17  0.21
luces_s[[1]][2000:2010,2000:2010] %>% table() # Sustraer una parte de la matriz
## .
## 0.0599999986588955 0.0700000002980232 0.0799999982118607 0.0900000035762787 
##                  3                  3                  4                  3 
##  0.100000001490116  0.109999999403954  0.119999997317791  0.129999995231628 
##                  4                 10                 10                 15 
##  0.140000000596046  0.150000005960464  0.159999996423721  0.170000001788139 
##                  8                 14                 12                  7 
##  0.180000007152557  0.189999997615814  0.200000002980232  0.209999993443489 
##                  6                  9                  3                  4 
##  0.219999998807907  0.239999994635582               0.25  0.280000001192093 
##                  2                  2                  1                  1
## puedo reproyectar un raster?
st_crs(luces_s)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
luces_new_crs = st_transform(luces_s,crs=4126)
luces_s[[1]][2000:2010,2000:2010] # no se alteran las geometrias
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24  0.28  0.14
##  [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16  0.25  0.14
##  [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10  0.17  0.12
##  [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12  0.14  0.12
##  [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15  0.13  0.11
##  [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15  0.10  0.13
##  [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14  0.13  0.13
##  [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13  0.19  0.19
##  [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09  0.13  0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11  0.11  0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18  0.17  0.21
luces_new_crs[[1]][2000:2010,2000:2010] # no se alteran las geometrias
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
##  [1,] 0.15 0.12 0.09 0.17 0.14 0.12 0.18 0.20 0.24  0.28  0.14
##  [2,] 0.15 0.12 0.07 0.16 0.14 0.13 0.12 0.16 0.16  0.25  0.14
##  [3,] 0.06 0.10 0.13 0.20 0.15 0.11 0.10 0.09 0.10  0.17  0.12
##  [4,] 0.07 0.13 0.18 0.16 0.11 0.20 0.21 0.08 0.12  0.14  0.12
##  [5,] 0.08 0.13 0.16 0.16 0.15 0.22 0.21 0.19 0.15  0.13  0.11
##  [6,] 0.12 0.13 0.16 0.15 0.13 0.18 0.19 0.17 0.15  0.10  0.13
##  [7,] 0.12 0.15 0.19 0.24 0.19 0.22 0.15 0.11 0.14  0.13  0.13
##  [8,] 0.16 0.11 0.08 0.18 0.21 0.19 0.11 0.15 0.13  0.19  0.19
##  [9,] 0.13 0.12 0.19 0.16 0.13 0.15 0.19 0.15 0.09  0.13  0.17
## [10,] 0.15 0.15 0.16 0.16 0.11 0.06 0.07 0.14 0.11  0.11  0.18
## [11,] 0.11 0.13 0.16 0.17 0.14 0.08 0.06 0.17 0.18  0.17  0.21
## plot data
plot(luces_s)

2.3 Filtrar usando la gemoetría

## download boundary
quilla <- getbb(place_name = "Barranquilla, Colombia", 
                featuretype = "boundary:administrative", 
                format_out = "sf_polygon") %>% .[1,]

leaflet() %>%
addTiles() %>%
addPolygons(data=quilla)
## cliping
l_quilla_1 = st_crop(x = luces_s , y = quilla) # crop luces de Colombia con polygono de Medellin

l_quilla_1[[1]] %>% t()
##       [,1] [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
##  [1,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##  [2,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##  [3,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##  [4,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##  [5,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##  [6,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##  [7,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##  [8,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##  [9,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [10,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [11,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [12,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [13,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [14,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [15,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [16,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [17,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [18,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [19,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [20,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [21,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [22,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [23,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [24,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [25,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA  6.03
## [26,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA  2.96  3.27
## [27,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA  2.14  2.26
## [28,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA  1.99
## [29,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA  1.95
## [30,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [31,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA  2.07
## [32,]   NA   NA    NA    NA    NA    NA  2.20  3.50    NA    NA  2.04  2.17
## [33,]   NA   NA    NA    NA    NA  1.65  1.94  1.93  1.88  2.04  2.21  2.26
## [34,]   NA   NA    NA    NA  1.78  2.46  2.97  2.20  2.08  3.82  5.71  5.30
## [35,]   NA   NA    NA    NA    NA  2.48  3.33  2.79  3.26  6.79  9.48  9.97
## [36,]   NA   NA    NA    NA    NA    NA  4.50  8.17 12.00 13.49 11.40  9.58
## [37,]   NA   NA    NA    NA 17.51 15.10 21.29 22.38 19.02 11.72  7.46  7.91
## [38,]   NA 1.84  9.53 15.32 25.18 18.67 24.94 22.62 14.00  4.76  3.73  3.57
## [39,]   NA 4.45 11.74    NA    NA  8.18    NA  8.98  6.40  3.54  3.17  3.00
## [40,]   NA   NA    NA    NA    NA    NA    NA    NA  3.81  3.43  2.76  2.53
## [41,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [42,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [43,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [44,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [45,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [46,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [47,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
## [48,]   NA   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
##       [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]  [,23]  [,24]
##  [1,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
##  [2,]    NA    NA    NA    NA  0.45    NA    NA    NA    NA    NA     NA     NA
##  [3,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
##  [4,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
##  [5,]    NA    NA    NA    NA    NA  0.64    NA    NA    NA    NA     NA     NA
##  [6,]    NA    NA    NA    NA    NA  0.56    NA    NA    NA    NA     NA     NA
##  [7,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
##  [8,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
##  [9,]    NA    NA    NA    NA    NA    NA  0.78    NA    NA    NA     NA     NA
## [10,]    NA    NA    NA    NA    NA    NA  0.86    NA    NA    NA     NA     NA
## [11,]    NA    NA    NA    NA    NA    NA  0.89    NA    NA    NA     NA     NA
## [12,]    NA    NA    NA    NA    NA    NA    NA  1.04    NA    NA     NA     NA
## [13,]    NA    NA    NA    NA    NA    NA  1.21  1.21  1.34    NA     NA     NA
## [14,]    NA    NA    NA    NA    NA  1.20    NA    NA  2.18    NA     NA     NA
## [15,]    NA    NA    NA    NA    NA    NA    NA    NA    NA  3.12     NA     NA
## [16,]    NA  1.60    NA    NA    NA    NA    NA    NA    NA  5.09  14.52  30.40
## [17,]  3.30  6.73    NA    NA    NA    NA    NA    NA    NA    NA  32.05  72.30
## [18,] 13.80 16.35    NA    NA    NA    NA    NA    NA    NA 42.66     NA  94.43
## [19,] 32.75 40.53 34.68    NA 12.67  6.47    NA    NA    NA    NA  95.34  76.51
## [20,] 25.76 44.97 46.67    NA    NA    NA    NA    NA    NA    NA  62.60  51.57
## [21,]    NA    NA    NA    NA    NA    NA    NA    NA    NA 15.24  21.01  27.16
## [22,]    NA    NA    NA    NA    NA    NA    NA    NA    NA 40.83  47.05  64.18
## [23,]    NA    NA 40.23    NA 27.55    NA    NA    NA 75.98 71.70  85.57 105.79
## [24,]    NA 13.75 14.97 15.46 11.35 18.27 29.89 53.61 82.99 93.33 101.92 111.54
## [25,]  5.50  4.56  4.78  4.01  4.50  5.21 14.28 42.38 76.66 85.32  91.35  99.35
## [26,]  3.04  2.84  3.02  3.31  3.98  4.93 14.24 34.13 62.03 71.77  79.13  87.62
## [27,]  2.28  2.48  2.90  3.18  4.69  8.84 23.25 28.93 35.22 45.27  62.79  76.90
## [28,]  2.11  2.35  2.58  3.02 11.18 32.95 71.85 63.93 44.63 33.76  42.14  58.06
## [29,]  2.07  2.20  2.53  3.29 12.92 43.91 82.88 75.04 43.40 22.66  35.22  49.98
## [30,]  2.11  2.26  2.52  2.90  6.60 25.38 28.53 29.68 28.91 42.55  45.05  52.81
## [31,]  2.25  2.42  2.76  2.81  4.35 11.00 16.23 17.50 42.77 62.58  64.82  55.21
## [32,]  2.31  2.49  3.30  4.51  8.08 20.60 36.12 36.38 49.91 70.47  69.25  54.57
## [33,]  2.48  2.82  4.23  7.03 13.51 24.19 39.68 49.92 55.94 69.37  62.20  45.75
## [34,]  5.34  6.25  7.87 10.68 21.66 23.54 28.73 51.42 59.80 58.90  48.16  43.53
## [35,] 12.41 17.60 15.63 22.35 33.51 28.47 29.17 47.31 63.17 52.35  46.63  45.12
## [36,] 15.48 23.54 34.28 43.02 56.55 37.83 52.31 66.26 65.72 53.72  47.85  49.67
## [37,] 11.54 18.39 37.05 68.12 61.43 46.70 58.65 70.88 66.73 56.34  52.74  53.29
## [38,]  4.42  7.21 17.60 36.16 54.00 38.88 30.27 26.63 40.38 57.27  55.93  56.22
## [39,]  3.14  3.63  9.59 16.15 18.85 13.85 10.31 12.22 26.62 54.40  56.93  57.73
## [40,]  2.69  3.13  4.20  7.02  6.92  9.13 12.85 18.41 26.70 60.24  69.57  64.83
## [41,]  2.67  5.43  8.80  8.01  9.52 17.38 25.22 37.87 40.71 54.39  69.89  67.61
## [42,]    NA    NA    NA 39.28    NA    NA 31.34 27.85 32.55 47.97  61.82  63.79
## [43,]    NA    NA    NA    NA    NA    NA 17.96  7.92 10.87 25.88  43.63  49.79
## [44,]    NA    NA    NA    NA    NA    NA    NA  6.99  8.95 30.74  47.43  54.56
## [45,]    NA    NA    NA    NA    NA    NA    NA    NA 10.48 28.05  47.61  58.19
## [46,]    NA    NA    NA    NA    NA    NA    NA    NA  6.33 13.28  27.71  47.00
## [47,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA  12.52  35.70
## [48,]    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA     NA     NA
##        [,25]  [,26] [,27]  [,28] [,29]  [,30] [,31] [,32] [,33] [,34]  [,35]
##  [1,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##  [2,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##  [3,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##  [4,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##  [5,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##  [6,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##  [7,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##  [8,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##  [9,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [10,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [11,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [12,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [13,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [14,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [15,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [16,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [17,]  44.35     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [18,]  92.54  72.77 38.77     NA    NA     NA    NA    NA    NA    NA     NA
## [19,]  56.01  59.52 47.77  39.21    NA     NA    NA    NA    NA    NA     NA
## [20,]  28.10  32.52 63.07  76.49 85.56  53.18    NA    NA    NA    NA     NA
## [21,]  46.55  50.28 71.32  83.99 85.30  63.33 52.52    NA    NA    NA     NA
## [22,]  76.44  77.28 78.88  65.12 60.87  54.69 53.39 43.20    NA    NA     NA
## [23,]  99.95  84.11 79.16  66.23 53.39  47.17 49.70 51.41 35.46 22.53     NA
## [24,] 102.47  76.40 81.48  79.32 66.51  56.38 59.51 62.07 46.97 25.86  10.23
## [25,]  96.30  91.31 86.57  86.37 84.80  69.04 69.33 66.04 51.06 27.09  19.61
## [26,]  95.27 100.51 87.12  86.48 85.51  77.50 65.75 62.99 51.63 51.20  36.54
## [27,]  89.58  96.36 92.16  90.84 87.09  82.91 75.94 66.18 67.95 74.07  70.44
## [28,]  68.38  83.42 87.67  95.82 90.66  76.26 78.83 78.96 75.55 81.83  70.24
## [29,]  60.20  68.46 79.86 100.32 91.54  73.72 75.44 93.24 92.60 77.20  65.92
## [30,]  58.82  61.30 65.44  69.87 71.60  71.90 75.94 88.63 98.86 86.89  71.64
## [31,]  55.67  55.18 55.81  61.82 66.07  71.64 73.51 78.55 92.67 85.76  82.93
## [32,]  48.39  55.60 57.99  59.12 65.07  69.80 67.73 73.24 83.59 81.23  83.38
## [33,]  43.38  48.67 58.43  58.08 60.75  62.08 58.91 66.09 72.25 70.98  76.85
## [34,]  42.76  48.42 56.79  54.78 54.65  60.34 60.78 67.11 64.15 64.04  64.97
## [35,]  50.09  56.28 57.06  53.92 57.21  72.60 79.05 76.00 68.02 64.23  63.44
## [36,]  57.15  60.50 53.60  50.88 54.49  66.49 70.52 66.77 71.46 68.73  68.48
## [37,]  67.56  65.18 53.31  55.93 60.69  63.25 60.49 60.49 71.85 89.37  73.30
## [38,]  63.81  57.03 52.76  59.40 61.73  75.28 80.32 70.60 82.59 84.18  70.54
## [39,]  57.26  51.44 49.65  50.65 51.92  69.87 86.02 78.78 92.53 87.69  74.86
## [40,]  56.72  47.53 44.60  44.62 49.16  59.19 73.68 80.36 92.06 99.48  96.92
## [41,]  56.60  48.08 45.28  47.42 49.30  67.65 71.78 69.61 84.80 96.47 111.94
## [42,]  62.72  59.50 49.69  52.51 58.76  80.24 76.63 67.45 72.62 79.86  97.51
## [43,]  63.71  60.43 52.18  59.41 72.81  93.18 94.09 87.39    NA    NA     NA
## [44,]  61.69  58.86 60.99  69.98 79.17  96.70 97.20    NA    NA    NA     NA
## [45,]  64.35  67.04 71.71  82.92 91.30 119.18    NA    NA    NA    NA     NA
## [46,]  57.68     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [47,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
## [48,]     NA     NA    NA     NA    NA     NA    NA    NA    NA    NA     NA
##       [,36] [,37] [,38] [,39] [,40] [,41]
##  [1,]    NA    NA    NA    NA    NA    NA
##  [2,]    NA    NA    NA    NA    NA    NA
##  [3,]    NA    NA    NA    NA    NA    NA
##  [4,]    NA    NA    NA    NA    NA    NA
##  [5,]    NA    NA    NA    NA    NA    NA
##  [6,]    NA    NA    NA    NA    NA    NA
##  [7,]    NA    NA    NA    NA    NA    NA
##  [8,]    NA    NA    NA    NA    NA    NA
##  [9,]    NA    NA    NA    NA    NA    NA
## [10,]    NA    NA    NA    NA    NA    NA
## [11,]    NA    NA    NA    NA    NA    NA
## [12,]    NA    NA    NA    NA    NA    NA
## [13,]    NA    NA    NA    NA    NA    NA
## [14,]    NA    NA    NA    NA    NA    NA
## [15,]    NA    NA    NA    NA    NA    NA
## [16,]    NA    NA    NA    NA    NA    NA
## [17,]    NA    NA    NA    NA    NA    NA
## [18,]    NA    NA    NA    NA    NA    NA
## [19,]    NA    NA    NA    NA    NA    NA
## [20,]    NA    NA    NA    NA    NA    NA
## [21,]    NA    NA    NA    NA    NA    NA
## [22,]    NA    NA    NA    NA    NA    NA
## [23,]    NA    NA    NA    NA    NA    NA
## [24,]    NA    NA    NA    NA    NA    NA
## [25,] 12.84    NA    NA    NA    NA    NA
## [26,] 31.83 16.36    NA    NA    NA    NA
## [27,] 46.73 21.85    NA    NA    NA    NA
## [28,] 36.41 18.83 11.97    NA    NA    NA
## [29,] 41.21 23.91 22.20    NA    NA    NA
## [30,] 60.35 49.21 42.73 29.10    NA    NA
## [31,] 82.86 69.00 61.91 51.15    NA    NA
## [32,] 95.06 82.41 64.15 52.33 45.85    NA
## [33,] 87.11 81.52 44.85 37.96 39.19    NA
## [34,] 68.56 65.19 51.40 72.00 66.22    NA
## [35,] 70.97 66.07 56.48 89.83 82.82    NA
## [36,] 75.41 76.00 55.55 51.73 47.74    NA
## [37,] 72.00 73.48 63.38 50.81 37.12    NA
## [38,] 71.74 76.27 80.48 83.79 80.64    NA
## [39,] 83.80 94.17 93.37 85.02 82.15    NA
## [40,] 87.09 91.57 88.35 76.77 37.66 21.85
## [41,]    NA    NA    NA    NA    NA    NA
## [42,]    NA    NA    NA    NA    NA    NA
## [43,]    NA    NA    NA    NA    NA    NA
## [44,]    NA    NA    NA    NA    NA    NA
## [45,]    NA    NA    NA    NA    NA    NA
## [46,]    NA    NA    NA    NA    NA    NA
## [47,]    NA    NA    NA    NA    NA    NA
## [48,]    NA    NA    NA    NA    NA    NA
## Plot data
ggplot() + geom_stars(data=l_quilla_1 , aes(y=y,x=x,fill=date_202301)) + # plot raster
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_bw() 

[3.] Aplicación: Actividad económica en Barranquilla

3.1 Importar datos 201301

## load data
l_quilla_0 = read_stars("input/raster/night_light_201301.tif") %>% st_crop(quilla)
names(l_quilla_0) = "date_201301"

ggplot() + geom_stars(data=l_quilla_0 , aes(y=y,x=x,fill=date_201301)) + # plot raster
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_bw() 

## stack rasters: adicionar bandas
l_quilla = c(l_quilla_0,l_quilla_1)
l_quilla
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##              Min. 1st Qu. Median     Mean 3rd Qu.   Max. NA's
## date_201301  0.15  6.1175 48.555 42.78252 69.3050 135.63 1270
## date_202301  0.45 15.2600 51.175 46.05638 69.8525 119.18 1270
## dimension(s):
##   from   to offset     delta refsys point x/y
## x  982 1022 -79.01  0.004167 WGS 84 FALSE [x]
## y  325  372  12.46 -0.004167 WGS 84 FALSE [y]

3.2 Raster a datos vectoriales

## st_as_sf
puntos_quilla = st_as_sf(x = l_quilla, as_points = T, na.rm = T) # raster to sf (points)
puntos_quilla
## Simple feature collection with 698 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.92083 ymin: 10.9125 xmax: -74.75417 ymax: 11.10833
## Geodetic CRS:  WGS 84
## First 10 features:
##    date_201301 date_202301                   geometry
## 1         0.15        0.45 POINT (-74.85417 11.10417)
## 2         0.26        0.64    POINT (-74.85 11.09167)
## 3         0.27        0.56     POINT (-74.85 11.0875)
## 4         0.32        0.78   POINT (-74.84583 11.075)
## 5         0.42        0.86 POINT (-74.84583 11.07083)
## 6         0.39        0.89 POINT (-74.84583 11.06667)
## 7         0.92        1.04  POINT (-74.84167 11.0625)
## 8         0.91        1.21 POINT (-74.84583 11.05833)
## 9         1.25        1.21 POINT (-74.84167 11.05833)
## 10        1.13        1.34  POINT (-74.8375 11.05833)
poly_quilla = st_as_sf(x = l_quilla, as_points = F, na.rm = T) # raster to sf (polygons)
poly_quilla
## Simple feature collection with 698 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.92292 ymin: 10.91042 xmax: -74.75208 ymax: 11.11042
## Geodetic CRS:  WGS 84
## First 10 features:
##    date_201301 date_202301                       geometry
## 1         0.15        0.45 POLYGON ((-74.85625 11.1062...
## 2         0.26        0.64 POLYGON ((-74.85208 11.0937...
## 3         0.27        0.56 POLYGON ((-74.85208 11.0895...
## 4         0.32        0.78 POLYGON ((-74.84792 11.0770...
## 5         0.42        0.86 POLYGON ((-74.84792 11.0729...
## 6         0.39        0.89 POLYGON ((-74.84792 11.0687...
## 7         0.92        1.04 POLYGON ((-74.84375 11.0645...
## 8         0.91        1.21 POLYGON ((-74.84792 11.0604...
## 9         1.25        1.21 POLYGON ((-74.84375 11.0604...
## 10        1.13        1.34 POLYGON ((-74.83958 11.0604...
## Plot data
ggplot() + 
geom_sf(data = puntos_quilla , aes(col=date_202301))  + 
scale_color_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_test() 

ggplot() + 
geom_sf(data = poly_quilla , aes(fill=date_202301),col=NA)  + 
scale_fill_viridis(option="A" , na.value='white') +
geom_sf(data=quilla , fill=NA , col="green") + theme_test()

3.3 Asignar el valor del pixel a un lugar

## get polygon
puerto <- getbb(place_name = "Puerto de Barranquilla - Sociedad Portuaria, Colombia", 
                featuretype = "barrier:wall", 
                format_out = "sf_polygon")

leaflet() %>%
addTiles() %>%
addPolygons(data=puerto)
## asignar luminosidad
l_puerto = st_join(puerto,poly_quilla)

## variacion promedio 
df = l_puerto
st_geometry(df) = NULL

df %>% 
summarise(pre=mean(date_201301,na.rm=T) , 
          post=mean(date_202301,na.rm=T)) %>%
mutate(ratio=post/pre-1)
##      pre   post     ratio
## 1 52.647 59.537 0.1308717

[4.] Aplicación: Cambio de cobertura de suelo?

## read raster
land_cover = c(read_stars("input/raster/discreta_2015.tif"),
               read_stars("input/raster/discreta_2019.tif"))
0.000992063*111000 # resolucion
## [1] 110.119
## read polygon
quilla_buffer <- st_buffer(x = quilla,dist = 2000)

leaflet() %>%
addTiles() %>%
addPolygons(data=quilla_buffer)
quilla_border <- st_difference(x = quilla_buffer , y = quilla)

leaflet() %>%
addTiles() %>%
addPolygons(data=quilla_border)
## cliping raster
quilla_cover = land_cover %>% st_crop(quilla_buffer)
plot(quilla_cover) ## view data

## cliping raster
border = land_cover %>% st_crop(quilla_border)
plot(border) ## view data

## raster to sf
border_sf = st_as_sf(x=border, as_points = F, na.rm = T) 

## cuantos pixeles se urbanizaron?
border_sf = border_sf %>% 
            rename(c_2015=discreta_2015.tif,c_2019=discreta_2019.tif)
border_sf = border_sf %>% 
            mutate(build = ifelse(c_2015!=50 & c_2019==50,1,0))
table(border_sf$build)
## 
##     0     1 
## 13990    46

[5.] Aplicación: Clusteres espaciales

## restaurantes
points <- opq(bbox = getbb("Bogota Colombia")) %>%
          add_osm_feature(key = "amenity", value = "restaurant") %>%
          osmdata_sf() %>% .$osm_points %>% select(osm_id,name) %>% mutate(restaurant=1)

leaflet() %>% addTiles() %>% addCircles(data=points)
## download boundary
city <- getbb(place_name = "Bogota, Colombia", 
        featuretype = "boundary:administrative", 
        format_out = "sf_polygon") %>% .[[2]]

leaflet() %>% addTiles() %>% addPolygons(data=city)
## afine transformations
points = st_transform(points,3857)
city = st_transform(city,3857)

## sf to ppp
points_p <- as.ppp(X = points["restaurant"])
Window(points_p) <- as.owin(city)
plot(points_p)

#--------------#
## summary
summary(points_p) # las unidades son en metros
## Marked planar point pattern:  6882 points
## Average intensity 1.733022e-05 points per square unit
## 
## Coordinates are given to 10 decimal places
## 
## marks are numeric, of type 'double'
## Summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1 
## 
## Window: polygonal boundary
## 5 separate polygons (no holes)
##            vertices        area relative.area
## polygon 1      7590 395161000.0      0.995000
## polygon 2       276   1524740.0      0.003840
## polygon 3        72    235244.0      0.000592
## polygon 4        45     75844.4      0.000191
## polygon 5        81    112582.0      0.000284
## enclosing rectangle: [-8262524, -8238784] x [498234.9, 538674.6] units
##                      (23740 x 40440 units)
## Window area = 397110000 square units
## Fraction of frame area: 0.414
intensity(points_p)
## [1] 1.733022e-05
## reescalar
points_p = rescale(X = points_p , s = 100 , unitname = "100-metros")
unitname(points_p)
## 100-metros / 100-metros
summary(points_p)
## Marked planar point pattern:  6882 points
## Average intensity 0.1733022 points per square 100-metros
## 
## Coordinates are given to 12 decimal places
## 
## marks are numeric, of type 'double'
## Summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1 
## 
## Window: polygonal boundary
## 5 separate polygons (no holes)
##            vertices        area relative.area
## polygon 1      7590 39516.10000      0.995000
## polygon 2       276   152.47400      0.003840
## polygon 3        72    23.52440      0.000592
## polygon 4        45     7.58444      0.000191
## polygon 5        81    11.25820      0.000284
## enclosing rectangle: [-82625.24, -82387.84] x [4982.349, 5386.746] 100-metros
##                      (237.4 x 404.4 100-metros)
## Window area = 39711 square 100-metros
## Unit of length: 1 100-metros
## Fraction of frame area: 0.414
intensity(points_p)
## [1] 0.1733022
## estimated standard error for intensity
sqrt(intensity(points_p)/summary(points_p)$window$area)
## [1] 0.00208904
## plot 
plot(density(points_p,5)) ## plot density

contour(density(points_p,bw.ppl(points_p)), axes = FALSE) ## contour density

#--------------#
## quadrants
points_q = quadratcount(points_p, nx = 4, ny = 6)
points_q
## tile
## Tile row 1, col 3 Tile row 1, col 4 Tile row 2, col 2 Tile row 2, col 3 
##                 0                20               100               345 
## Tile row 2, col 4 Tile row 3, col 1 Tile row 3, col 2 Tile row 3, col 3 
##               739                 0               373               861 
## Tile row 3, col 4 Tile row 4, col 1 Tile row 4, col 2 Tile row 4, col 3 
##               864                53               647              2029 
## Tile row 4, col 4 Tile row 5, col 1 Tile row 5, col 2 Tile row 5, col 3 
##               263                17               183               377 
## Tile row 5, col 4 Tile row 6, col 2 Tile row 6, col 3 
##                 0                 0                11
plot(points_q, cex = 2 , col="red")

Referencias

  • Lovelace, R., Nowosad, J., & Muenchow, J. (2019). Geocomputation with R. [Ver aquí]

    • Cap. 4: Spatial Data Operations
    • Cap. 5: Geometry Operations
    • Cap. 6: Reprojecting geographic data
    • Cap. 11: Statistical learning
  • Bivand, R. S., Pebesma, E. J., Gómez-Rubio, V., & Pebesma, E. J. (2013). Applied spatial data analysis with R. [Ver aquí]

    • Cap. 4: Spatial Data Import and Export
    • Cap. 7: Spatial Point Pattern Analysis
    • Cap. 8: Interpolation and Geostatistics